{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f77bc74d-58a1-4c14-b477-f52d28f2a869",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| default_exp models.mlp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15392f6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "12fa25a4",
   "metadata": {},
   "source": [
    "# MLP\n",
    "> One of the simplest neural architectures are Multi Layer Perceptrons (`MLP`) composed of stacked Fully Connected Neural Networks trained with backpropagation. Each node in the architecture is capable of modeling non-linear relationships granted by their activation functions. Novel activations like Rectified Linear Units (`ReLU`) have greatly improved the ability to fit deeper networks overcoming gradient vanishing problems that were associated with `Sigmoid` and `TanH` activations. For the forecasting task the last layer is changed to follow a auto-regression problem.<br><br>**References**<br>-[Rosenblatt, F. (1958). \"The perceptron: A probabilistic model for information storage and organization in the brain.\"](https://psycnet.apa.org/record/1959-09865-001)<br>-[Fukushima, K. (1975). \"Cognitron: A self-organizing multilayered neural network.\"](https://pascal-francis.inist.fr/vibad/index.php?action=getRecordDetail&idt=PASCAL7750396723)<br>-[Vinod Nair, Geoffrey E. Hinton (2010). \"Rectified Linear Units Improve Restricted Boltzmann Machines\"](https://www.cs.toronto.edu/~fritz/absps/reluICML.pdf)<br>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e6036ce9",
   "metadata": {},
   "source": [
    "![Figure 1. Three layer MLP with autorregresive inputs.](imgs_models/mlp.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2508f7a9-1433-4ad8-8f2f-0078c6ed6c3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "from fastcore.test import test_eq\n",
    "from nbdev.showdoc import show_doc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44065066-e72a-431f-938f-1528adef9fe8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "from typing import Optional\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "\n",
    "from neuralforecast.losses.pytorch import MAE\n",
    "from neuralforecast.common._base_windows import BaseWindows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce70cd14-ecb1-4205-8511-fecbd26c8408",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "class MLP(BaseWindows):\n",
    "    \"\"\" MLP\n",
    "\n",
    "    Simple Multi Layer Perceptron architecture (MLP). \n",
    "    This deep neural network has constant units through its layers, each with\n",
    "    ReLU non-linearities, it is trained using ADAM stochastic gradient descent.\n",
    "    The network accepts static, historic and future exogenous data, flattens \n",
    "    the inputs and learns fully connected relationships against the target variable.\n",
    "\n",
    "    **Parameters:**<br>\n",
    "    `h`: int, forecast horizon.<br>\n",
    "    `input_size`: int, considered autorregresive inputs (lags), y=[1,2,3,4] input_size=2 -> lags=[1,2].<br>\n",
    "    `stat_exog_list`: str list, static exogenous columns.<br>\n",
    "    `hist_exog_list`: str list, historic exogenous columns.<br>\n",
    "    `futr_exog_list`: str list, future exogenous columns.<br>\n",
    "    `exclude_insample_y`: bool=False, the model skips the autoregressive features y[t-input_size:t] if True.<br>\n",
    "    `n_layers`: int, number of layers for the MLP.<br>\n",
    "    `hidden_size`: int, number of units for each layer of the MLP.<br>\n",
    "    `loss`: PyTorch module, instantiated train loss class from [losses collection](https://nixtla.github.io/neuralforecast/losses.pytorch.html).<br>\n",
    "    `valid_loss`: PyTorch module=`loss`, instantiated valid loss class from [losses collection](https://nixtla.github.io/neuralforecast/losses.pytorch.html).<br>\n",
    "    `max_steps`: int=1000, maximum number of training steps.<br>\n",
    "    `learning_rate`: float=1e-3, Learning rate between (0, 1).<br>\n",
    "    `num_lr_decays`: int=-1, Number of learning rate decays, evenly distributed across max_steps.<br>\n",
    "    `early_stop_patience_steps`: int=-1, Number of validation iterations before early stopping.<br>\n",
    "    `val_check_steps`: int=100, Number of training steps between every validation loss check.<br>\n",
    "    `batch_size`: int=32, number of different series in each batch.<br>\n",
    "    `valid_batch_size`: int=None, number of different series in each validation and test batch, if None uses batch_size.<br>\n",
    "    `windows_batch_size`: int=1024, number of windows to sample in each training batch, default uses all.<br>\n",
    "    `inference_windows_batch_size`: int=-1, number of windows to sample in each inference batch, -1 uses all.<br>\n",
    "    `start_padding_enabled`: bool=False, if True, the model will pad the time series with zeros at the beginning, by input size.<br>\n",
    "    `step_size`: int=1, step size between each window of temporal data.<br>\n",
    "    `scaler_type`: str='identity', type of scaler for temporal inputs normalization see [temporal scalers](https://nixtla.github.io/neuralforecast/common.scalers.html).<br>\n",
    "    `random_seed`: int=1, random_seed for pytorch initializer and numpy generators.<br>\n",
    "    `num_workers_loader`: int=os.cpu_count(), workers to be used by `TimeSeriesDataLoader`.<br>\n",
    "    `drop_last_loader`: bool=False, if True `TimeSeriesDataLoader` drops last non-full batch.<br>\n",
    "    `alias`: str, optional,  Custom name of the model.<br>\n",
    "    `**trainer_kwargs`: int,  keyword trainer arguments inherited from [PyTorch Lighning's trainer](https://pytorch-lightning.readthedocs.io/en/stable/api/pytorch_lightning.trainer.trainer.Trainer.html?highlight=trainer).<br>    \n",
    "    \"\"\"\n",
    "    # Class attributes\n",
    "    SAMPLING_TYPE = 'windows'\n",
    "    \n",
    "    def __init__(self,\n",
    "                 h,\n",
    "                 input_size,\n",
    "                 futr_exog_list = None,\n",
    "                 hist_exog_list = None,\n",
    "                 stat_exog_list = None,\n",
    "                 exclude_insample_y = False,\n",
    "                 num_layers = 2,\n",
    "                 hidden_size = 1024,\n",
    "                 loss = MAE(),\n",
    "                 valid_loss = None,\n",
    "                 max_steps: int = 1000,\n",
    "                 learning_rate: float = 1e-3,\n",
    "                 num_lr_decays: int = -1,\n",
    "                 early_stop_patience_steps: int =-1,\n",
    "                 val_check_steps: int = 100,\n",
    "                 batch_size: int = 32,\n",
    "                 valid_batch_size: Optional[int] = None,\n",
    "                 windows_batch_size = 1024,\n",
    "                 inference_windows_batch_size = -1,\n",
    "                 start_padding_enabled = False,\n",
    "                 step_size: int = 1,\n",
    "                 scaler_type: str = 'identity',\n",
    "                 random_seed: int = 1,\n",
    "                 num_workers_loader: int = 0,\n",
    "                 drop_last_loader: bool = False,\n",
    "                 **trainer_kwargs):\n",
    "\n",
    "        # Inherit BaseWindows class\n",
    "        super(MLP, self).__init__(h=h,\n",
    "                                  input_size=input_size,\n",
    "                                  futr_exog_list=futr_exog_list,\n",
    "                                  hist_exog_list=hist_exog_list,\n",
    "                                  stat_exog_list=stat_exog_list,\n",
    "                                  exclude_insample_y = exclude_insample_y,\n",
    "                                  loss=loss,\n",
    "                                  valid_loss=valid_loss,\n",
    "                                  max_steps=max_steps,\n",
    "                                  learning_rate=learning_rate,\n",
    "                                  num_lr_decays=num_lr_decays,\n",
    "                                  early_stop_patience_steps=early_stop_patience_steps,\n",
    "                                  val_check_steps=val_check_steps,\n",
    "                                  batch_size=batch_size,\n",
    "                                  valid_batch_size=valid_batch_size,\n",
    "                                  windows_batch_size=windows_batch_size,\n",
    "                                  inference_windows_batch_size=inference_windows_batch_size,\n",
    "                                  start_padding_enabled=start_padding_enabled,\n",
    "                                  step_size=step_size,\n",
    "                                  scaler_type=scaler_type,\n",
    "                                  num_workers_loader=num_workers_loader,\n",
    "                                  drop_last_loader=drop_last_loader,\n",
    "                                  random_seed=random_seed,\n",
    "                                  **trainer_kwargs)\n",
    "\n",
    "        # Architecture\n",
    "        self.num_layers = num_layers\n",
    "        self.hidden_size = hidden_size\n",
    "\n",
    "        self.futr_input_size = len(self.futr_exog_list)\n",
    "        self.hist_input_size = len(self.hist_exog_list)\n",
    "        self.stat_input_size = len(self.stat_exog_list)\n",
    "\n",
    "        input_size_first_layer = input_size + self.hist_input_size * input_size + \\\n",
    "                                 self.futr_input_size*(input_size + h) + self.stat_input_size\n",
    "\n",
    "        # MultiLayer Perceptron\n",
    "        layers = [nn.Linear(in_features=input_size_first_layer, out_features=hidden_size)]\n",
    "        for i in range(num_layers - 1):\n",
    "            layers += [nn.Linear(in_features=hidden_size, out_features=hidden_size)]\n",
    "        self.mlp = nn.ModuleList(layers)\n",
    "\n",
    "        # Adapter with Loss dependent dimensions\n",
    "        self.out = nn.Linear(in_features=hidden_size, \n",
    "                             out_features=h * self.loss.outputsize_multiplier)\n",
    "\n",
    "    def forward(self, windows_batch):\n",
    "\n",
    "        # Parse windows_batch\n",
    "        insample_y    = windows_batch['insample_y']\n",
    "        futr_exog     = windows_batch['futr_exog']\n",
    "        hist_exog     = windows_batch['hist_exog']\n",
    "        stat_exog     = windows_batch['stat_exog']\n",
    "\n",
    "        # Flatten MLP inputs [B, L+H, C] -> [B, (L+H)*C]\n",
    "        # Contatenate [ Y_t, | X_{t-L},..., X_{t} | F_{t-L},..., F_{t+H} | S ]\n",
    "        batch_size = len(insample_y)\n",
    "        if self.hist_input_size > 0:\n",
    "            insample_y = torch.cat(( insample_y, hist_exog.reshape(batch_size,-1) ), dim=1)\n",
    "\n",
    "        if self.futr_input_size > 0:\n",
    "            insample_y = torch.cat(( insample_y, futr_exog.reshape(batch_size,-1) ), dim=1)\n",
    "\n",
    "        if self.stat_input_size > 0:\n",
    "            insample_y = torch.cat(( insample_y, stat_exog.reshape(batch_size,-1) ), dim=1)\n",
    "\n",
    "        y_pred = insample_y.clone()\n",
    "        for layer in self.mlp:\n",
    "             y_pred = torch.relu(layer(y_pred))\n",
    "        y_pred = self.out(y_pred)\n",
    "\n",
    "        y_pred = y_pred.reshape(batch_size, self.h, \n",
    "                                self.loss.outputsize_multiplier)\n",
    "        y_pred = self.loss.domain_map(y_pred)\n",
    "        return y_pred"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfc06a06",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(MLP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a23696b",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(MLP.fit, name='MLP.fit')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8475d33",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(MLP.predict, name='MLP.predict')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac34472d-5670-45b5-a7c5-1dba54f8e782",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "import logging\n",
    "import warnings\n",
    "logging.getLogger(\"pytorch_lightning\").setLevel(logging.ERROR)\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8b06414f-117e-42ce-abad-46ad4b8372b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "# test performance fit/predict method\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from neuralforecast.utils import AirPassengersDF as Y_df\n",
    "from neuralforecast.tsdataset import TimeSeriesDataset\n",
    "\n",
    "\n",
    "Y_train_df = Y_df[Y_df.ds<='1959-12-31'] # 132 train\n",
    "Y_test_df = Y_df[Y_df.ds>'1959-12-31']   # 12 test\n",
    "\n",
    "dataset, *_ = TimeSeriesDataset.from_df(Y_train_df)\n",
    "model = MLP(h=12, input_size=24, max_steps=1)\n",
    "model.fit(dataset=dataset)\n",
    "y_hat = model.predict(dataset=dataset)\n",
    "Y_test_df['MLP'] = y_hat\n",
    "\n",
    "#test we recover the same forecast\n",
    "y_hat2 = model.predict(dataset=dataset)\n",
    "test_eq(y_hat, y_hat2)\n",
    "\n",
    "pd.concat([Y_train_df, Y_test_df]).drop('unique_id', axis=1).set_index('ds').plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30f70ba2-2c57-4710-808b-727f7e0169a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "# test no leakage with test_size\n",
    "dataset, *_ = TimeSeriesDataset.from_df(Y_df)\n",
    "model = MLP(h=12, input_size=24, max_steps=1)\n",
    "model.fit(dataset=dataset, test_size=12)\n",
    "y_hat_test = model.predict(dataset=dataset, step_size=1)\n",
    "np.testing.assert_almost_equal(\n",
    "    y_hat, \n",
    "    y_hat_test,\n",
    "    decimal=4\n",
    ")\n",
    "# test we recover the same forecast\n",
    "y_hat_test2 = model.predict(dataset=dataset, step_size=1)\n",
    "test_eq(y_hat_test, y_hat_test2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3a481bb-65df-444f-9ade-ec23c71d5304",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "# test validation step\n",
    "dataset, *_ = TimeSeriesDataset.from_df(Y_train_df)\n",
    "model = MLP(h=12, input_size=24, step_size=1, \n",
    "            hidden_size=1024, num_layers=2,\n",
    "            max_steps=1)\n",
    "model.fit(dataset=dataset, val_size=12)\n",
    "y_hat_w_val = model.predict(dataset=dataset)\n",
    "Y_test_df['MLP'] = y_hat_w_val\n",
    "\n",
    "pd.concat([Y_train_df, Y_test_df]).drop('unique_id', axis=1).set_index('ds').plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51d5df26-2bc7-405b-a421-2cace4a39d10",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "# test no leakage with test_size and val_size\n",
    "dataset, *_ = TimeSeriesDataset.from_df(Y_df)\n",
    "model = MLP(h=12, input_size=24, step_size=1, \n",
    "            hidden_size=1024, num_layers=2,\n",
    "            max_steps=1)\n",
    "model.fit(dataset=dataset, val_size=12, test_size=12)\n",
    "y_hat_test_w_val = model.predict(dataset=dataset, step_size=1)\n",
    "np.testing.assert_almost_equal(y_hat_test_w_val,\n",
    "                               y_hat_w_val, decimal=4)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "9c61645f",
   "metadata": {},
   "source": [
    "## Usage Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72b60ba0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| eval: false\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pytorch_lightning as pl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from neuralforecast import NeuralForecast\n",
    "from neuralforecast.models import MLP\n",
    "from neuralforecast.losses.pytorch import MQLoss, DistributionLoss\n",
    "from neuralforecast.tsdataset import TimeSeriesDataset\n",
    "from neuralforecast.utils import AirPassengers, AirPassengersPanel, AirPassengersStatic\n",
    "\n",
    "Y_train_df = AirPassengersPanel[AirPassengersPanel.ds<AirPassengersPanel['ds'].values[-12]] # 132 train\n",
    "Y_test_df = AirPassengersPanel[AirPassengersPanel.ds>=AirPassengersPanel['ds'].values[-12]].reset_index(drop=True) # 12 test\n",
    "\n",
    "model = MLP(h=12, input_size=24,\n",
    "            loss=DistributionLoss(distribution='Normal', level=[80, 90]),\n",
    "            scaler_type='robust',\n",
    "            learning_rate=1e-3,\n",
    "            max_steps=200,\n",
    "            val_check_steps=10,\n",
    "            early_stop_patience_steps=2)\n",
    "\n",
    "fcst = NeuralForecast(\n",
    "    models=[model],\n",
    "    freq='M'\n",
    ")\n",
    "fcst.fit(df=Y_train_df, static_df=AirPassengersStatic, val_size=12)\n",
    "forecasts = fcst.predict(futr_df=Y_test_df)\n",
    "\n",
    "Y_hat_df = forecasts.reset_index(drop=False).drop(columns=['unique_id','ds'])\n",
    "plot_df = pd.concat([Y_test_df, Y_hat_df], axis=1)\n",
    "plot_df = pd.concat([Y_train_df, plot_df])\n",
    "\n",
    "plot_df = plot_df[plot_df.unique_id=='Airline1'].drop('unique_id', axis=1)\n",
    "plt.plot(plot_df['ds'], plot_df['y'], c='black', label='True')\n",
    "plt.plot(plot_df['ds'], plot_df['MLP-median'], c='blue', label='median')\n",
    "plt.fill_between(x=plot_df['ds'][-12:], \n",
    "                 y1=plot_df['MLP-lo-90'][-12:].values, \n",
    "                 y2=plot_df['MLP-hi-90'][-12:].values,\n",
    "                 alpha=0.4, label='level 90')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.plot()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "python3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
